9 Abstract 



A system is described which allows buyers to define their preferences and sellers to define their capabilities, 
then determines which trading points maximize the utility of the buyer. The system suggests trades by 
exploiting the flexibilities and tradeoffs encoded by both parties, thus providing win-win trades. A second 
level of optimization ranks the trades with all suppliers, allowing the buyer to rapidly determine the best 
alternatives. The system allows for rich negotiation spaces and supports continuous, discrete, and range or 
interval decision factors. 
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A Distance Function for Range Variables 



The distance between buyer range variable n = (r t ^i) an d seller range variable r 3 — (rj,r 3 ) is related to 
the overlap defined as 

overlap (r^Tj) = / dx N(x\ r t )M{x\ r 3 ) 

where Af(x;r z ) is a Gaussian in x centered at = (r t + r % )/2 with standard deviation a z — a(fi — r t ) and 
a is a tunable parameter. The range of integration is bounded by the seller's capabilities. Evaluating the 
integral, we find 



over lap (r t , r 3 ) — 



2^2ntf +<rj) 



exp 



erf 



erf -\-a 3 fii - <7 3 fi 3 +af 



- erf 



-a] + (TjfH - a 3 n 3 - o\ 



By adjusting a we can control the hardness of the boundary at the ends of the range (smaller a corresponds to 
sharper dropoff outside the range). The maximum overlap will occur when fi 3 = fii and a 3 = 0. Substituting 
those values in the above equation, we get 

r ; rf ^ 

The distance d(r t ,r 3 ) is related to the overlap so that it varies from 0 (for identical range variables) to 
oo (for maximally distant range variables). Consequently, we define 



R(r lz r 3 ) - - In 



overlap (r t , r 3 ) 
maxOverlap 



B From Rankings to Numerical Values 

The numerical values in the tables Z(x t: x 3 ) can be formed from the buyer's rankings in a number of ways. 
Since we are minimizing distance, we want preferred options to contribute little Z and progressively less 
favorable options to have progressively larger Z values. To approach this, we first consider the problem of 
defining a Z* value normalized to lie between 0 and 1 and then take Z = — ln[l — Z*]. 

One possibility for Z f is a linear scaling. If d z is the integer indicating the number of values in V iy 
r\°^ {x 3 ) indicates the buyer-specified rank (smaller rank means more preferable) for the option Xi — a £ T>i 
conditional on x 3 , then we might set Z\x % = a,x 3 ) = {r\ a \v 3 ) - l)/(d t - 1). Alternatively, we might use 
a non-linear scaling and set Z{xi = a,v 3 ) = (1 - /? r ^ > ^^~ 1 )/(l - P d ^~ l ) for some P G [0, 1] (exponential) 
or set Z{v % — a,v 3 ) = ((r[ a \v 3 ) - l)/(R t — l)) n for some n > 0 (algebraic). The corresponding scaling 
functions for Z are given in Table 4. 
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type 


Z' 


linear 
exponential 
algebraic 


In 
-In 


In r-^ 

l-f3R,-i "I 


pT< a \v,)-l _ fiR,-l 



Table 4: Possible numerical ranking functions. 



For those allowed values of ^ which are left unranked by the buyer, we can assign them a common 
numerical value greater than the lowest ranked value. For example if r % (x 3 ) is the largest ranked value we 
might assign all remaining ^ values to have Z(^,^-) = r»(xj)/(d f - 1) in the linear case or Z{yc %1 >c i ) = 
(l-p ri ^)/(l~P d *- 1 ) or Z(^xj) = (1- l/{r % {x,) + 1))/(1 - l/d t ) in the non-linear cases. 



C Determination of Scaling Factors 

In this section we supply the details necessary to determine the scaling factors Q r and Q d . In the equations 
defining the average distances we have need of the integrals f v du exp[-(u - u)* C" 1 (u - u)] and J v du (u - 
u)*C _1 (u-u) exp[-(u-u)*C" 1 (u-u)]. In general, these integrations are not analytically tractable; however 
we will assume that the region of integration defined by x and x extends at least 2 standard deviations out 
in all directions so that the volume of integration is effectively infinite. In this case we can evaluate the 
integrals in closed form: 



X 



du exp[-(u - u)*C 1 (u - u)] & ^J- 



detC" 1 
and 



du (u - u/C-^u - u) exp[-(u - u^C-^u - u)] 



h 2 y detC' 1 ' 

To this point we have not been specific about what is meant by the average over the range variables. In 
fact, we will assume that, however we define this average, it is independent of what the buyer requested and 
thus can be removed from inside the sum over the discrete variables. In this case we will concentrate on 
setting the discrete and continuous components equal. We have (d c )/n c = {d d )/n d where 

(d d ) _ Qd Ev d d eM-Qdd d } J v du exp[-d c ] (d c ) _ 1 E v exp[-Q d d d ] f v du d c exp[-d c ] 

nd n d £ v exp[-Q d d d ]/ K du exp[-d c ] n c ~ n c £ v exp[-Q d d d ] f v du exp[-d c ] 



and cancelling denominators Q d is determined by solving 

■^d d exp[--Q d dd] / du exp[-cZ c ] = — J^exp[~Q d d d ] f dud c ex.p[-d c 
Jv n c ^ Jv 



Qd 
n d 
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This can be simplified further given our assumptions about the Gaussian integral: 

Qd d d exp[-Q d d d ] exp[-Q d d d ] 
n d 4" Vdet C-Hv) " 2 ^ ^det C-i(v) ' 

In most cases there will also be at least one finite root, which is what we need. If there are a multiplicity 
of roots, we select that root having the smallest value to ensure that the distribution denned by the utility 
function is as flat as possible. If there are no real roots, then we select the Q d which minimizes the error. 

D Phase 1 Optimization 

In this section we detail the solution to the phase 1 optimization problem defined in Eq. (5). 

D.l General Formulation 

We begin by writing down the objective function: 

di(x,x) - (x<*> -M^'C" 1 ^ -Pi*)) 

where (see Eq. (8)) 

x« = (I-M(*)(^\{^ (8) 

as long as x { k s) (b\ 3 J) < x[ s) < x[ s) (b[ s J + 1). This must be optimized over x subject to the constraints 
Gi(V)x (s ) = gi(V) and G 2 0<)x< s ) = g 2 (*r) for a given jwt. Since >c is fixed throughout the optimization 
R(r; k) is just a constant which can be evaluated and then ignored as far as determination of x is concerned. 
We simplify our notation so that di(u) = (x< 5 ) - itfC" 1 (x^ - p) + ii(r). Similarly, 

X W = (l-M^({^ 

where we have defined the n c -vector s = (I - M^) _1 c and the n c x n' c matrix S = (I - M^) _1 M^\ 16 
n' c < n c is the number of continuous dimensions upon which the supplier response depends. Putting these 
pieces together we have 

dx(xW) = x< 6 >*Q({&W, {6{J}) X W + q'x^ + const (9) 

where Q is an n l c x n l c matrix given by 

Q({»i,l{ff2» = s*({&lt {C2})c- 1 s({ 6 « ! {&$}) 

16 In the next section we consider S further, taking advantage of any redundant variables to simplify the QP task. 
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q is an n'-vector given by 

and the constant term is 

const = (*({&«{&$}) - -tfCZ^sttbQM*}) -p) +R(r) 

It is easily shown that the matrix Q is positive semidefinite. di(u^) must be optimized subject to the 
constraints 

GiS({6W, {&g>»x<»> = gl - G lS ({6W, {&«}), (10) 
G 2 S({6ll < g 2 - G 3 s({6W, (11) 

x W « < s({6, ( :2, {&$}) + S({b[% {6<S})xW < x«(0 V t, (12) 
xW(i) <x<" < Vi (13) 

where xW(i)* = [^(fiW) ... and x(»>(i)' = [ Xl{b % + 1) ... ^(6^ + 1)] are the vec- 

tors of breakpoints for the f$ (k e [l,n c ]) functions. Similarly, for the buyer variables, x (6) (*)' = 
[xi(b^l) ■■■ z n {bfy\ and x< 6 >(i)' = [a?i (6^ + 1) ••• a; n (6^ + 1)] are the vectors of breakpoints 
for the /f^ (A; € [l,n c ]) functions. 

di is optimized in two steps. First we assume a value for the integer variables {6^} and fixing the 

constraints in Eqs. (10) and (11). We then minimize the quadratic 4 subject to these constraints using a 
standard quadratic programming (QP) algorithm 17 to find x<p\ ({&$}, {&$})• This result is then used to 
find x* (x; {b[% = s ({&«}, {&$}) + S({b[% {ig})^ ({6$}, {6$}) . By searching exhaustively 

over {b[*l} and {6^2} we determine 

x<*>(*)= argmin di (x*(x; {^2})^) 

Before turning to the manner in which this optimization is accomplished in section D.3, we discuss a few 
improvements to the input into the QP algorithm. 

D.2 The Quadratic Program 

In this section we discuss improvements in the specification of the quadratic programming task to make it 
better conditioned and more efficient. 

It is possible that the matrix Q used as input to the QP is singular or near singular, depending on the 
rank of Si. Consequently, it is important to use an interior point QP algorithm for which this is not a 
problem. 

17 Numerous third parties supply quadratic programming packages, some based on Interior-point methods which are appro- 
priate, in the present context where matrix Q may be singular. 
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D.2.1 Variable Elimination 



There may be rows or columns of S which are entirely zero. It is important to take advantage of this structure 
to eliminate unnecessary computation. 

If row % is zero, this indicates that supplier variable u[ s) is independent of any buyer variable and is simply 
equal to s % . If there are any constraints involving only supplier variable i, then these can be checked before 
calling the QP solver to eliminate wasteful computation. 18 If the constraint is satisfied, we can eliminate it 
and then call the QP function. If column j of S is all zero, this means that there is no dependence of u( s ) 
on u { p and so there is no need to optimize over that buyer variable. We will, however, still have to check 
that the constraints on are satisfied. This observation allows us to simplify the quadratic programming 
problem by eliminating this unnecessary variable. 

We formulate the problem to take advantage of both these observations. Assume there are m < n c 
non-zero rows of S and m } < n f c non-zero columns. Relabel the rows and columns so that 



S = 



Si 
0 



where Si is of size m x m! and, with a slight abuse of notation, the O's indicate matrices of the appropriate 
size. In like fashion we write 



4" 



x<*> = 



where x[ s) is an m vector, x£° is a,n n c - m vector, x^° is an ml vector, and x ( 2 b) is an n' c - m! vector. 
x^ s) indicates those components of x<» which vary with buyer variables while x^° indicates those variables 
which are constant. x< 6) are those buyer variables which affect at least one supplier variable and x^ 6) are 
those buyer variables which do not affect any supplier variables. If we also partition s as s* = [s* s |] then 

x£ tf) = Sl +Six[ 6) and x£° = s 2 . 

Finally, we will also have need of the partitioned covariance matrix 



^2,1 



^2,2 



where (C 1 1 2 . w ) t = C 2j {. w since C" 1 is symmetric- 
Including the dependencies on fi, the quadratic form we must optimize is 



WiCrjWx WiC^SWa 
WsC-jWi w 2 c 2 Jw 2 



si - /*i + Smf } 

S2 - M 2 



Additionally, the QP doesn't like to have a row of zeros in the constraint matrix. 
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where we have made the obvious definitions 

[Wi 0 

w = 

0 w 2 

Expanding the quadratic, the constant term is 



and ix = 



Mi 
M 2 



(si - AO'WxC! jW x ( Sl - Ml ) + (s 2 - M 2 ) f W 2 C 2 -jW 2 (s 2 - M2 ) + 2( Sl - Mi) t W 1 C^W 2 (s 2 - 
The linear term is 

q = ((si - Mi^WiCrjWiSj + 2(s 2 - M 2 ) i W 2 C 2 -}W 1 S 1 )x^ 
and the quadratic term is 

Q = (xf ) ) t S' 1 W 1 C-iW 1 S 1 xf ) . 

The constraints involving the optimization variables are 

G M Six^ = gi - Gi.isi - G 1;2 s 2 , 
G 2; iSixf ) < g 2 - G 2) isi - G 2;2 s 2 , 
d s) (i)<s 1 +S 1 4 b) <*{ s) (i) Vi, 
x^(i) <xf» <x< 6) (i) Vi 

where we have simplified notation by dropping the dependence on breakpoints and defined 

G, = [c t;1 G, ;2 ] 

i.e., G i;] is the columns of Gi involving the active buyer variables while G t;2 is formed from the inactive 
buyer variables. Similarly, the constraints which can be checked prior to even solving the QP are 



X2°W < s 2 < *2 } {i) Vi 



D.2.2 QP Normalization 



There are two types of normalizations we perform for increased numerical stability. If the scales of the 
continuous dimensions vary greatly (i.e., one dimension has an ideal value of 10~ 2 and another dimension 
has ideal value 10 6 ), then the matrix C _1 becomes badly conditioned, which makes S badly conditioned and 
thus makes the QP unstable. Similarly, if the breakpoints have radically, different lengths, the QP problem 
can also become unstable. We consider how to remedy each of these problems. 

Dimension normalization 
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The quadratic form to be minimized is (x - ^) t C- 1 (x - p). We modify this problem by defining the 
dimensionless variables through x t = in terms of which 

(x - ^'c-^x - n) = {<f> - lyc 1 ^ - 1) 

where C -1 = DC^D and D = diag(p). When making this variable transformation we must assure that 
the constraints on x are adjusted so that G 2 x = gl and G 2 x < g 2 become the constraints Q x 4> = gl and 
Gi$ < g 2 on <p where Q h2 = G 12 D. If there are global bounds of the form x < x < x, then these are 
translated to < 0 < 0 where = D _1 u and ~$ - D _I u. 

In addition to these changes, we also need to modify the capabilities appropriately. Recall Eq. (2) which 
reads (after dropping the dependence on the discrete variables) 

tf-E&W+fS^). d4) 

k 

In terms of the dimensionless variables we have 

~V - + - • (15) 

Thus, when we use the dimensionless variables, the function f hk needs to be modified so that the y dimension 
is divided by \x % and the x dimension is multiplied by fj, k . 

Once the optimization has determined an optimal 0, we return to the original variables through x = &<fi. 
Breakpoint range normalization 

For larger problems it is sometimes important to normalize the input variables for the QP solver. This is 
carried out in order to improve the QP search performance in situations where there is large variation in 
the range of the segments <p ( k b) (b k ) < ^ < 0< 6 >(& fe + 1) and the magnitude of their values. This can be 
accomplished by introducing the normalization variables <p {b} = (<f) k b) - 4> { k b) (b k ))/((f> k b) (b k + 1) - <^ 6) (&*))• 
Here will have minimum value of zero and maximum of one. We can now express 0< 6 > as 

4>W = A0 (6) +p (16) 
where A is the diagonal matrix A = diag(0< 6l+1) - <$ 1 \<$* +l) - ^ a) , ■ • • , ^£ m+1) - and p* - 

Expanding Eq.(9) using the definition of <j> in Eq.(16), we obtain 

di(w; {b hk }) -i0 (6) *ADQ({fe^})DA^ (6) + q'({& J>fc })DA^ + p*DQ({fe,, Jfe })DA^ ) + 
lp t DQ({&,,,})D P + q*({6 l>fc })Dp. 
Hence the quadratic programming problem to be solved in terms of <f> is 

QP(x«; {bi, k }) = ^tQdb^})^ + 4*({6 tii })u(*) > (i 7 ) 
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where Q = ADQ({6,, 4 })DA and q* = (q t ({&,,*})+p t DQ({&, i *}))DA. Similarly, the equality and inequality 
constraints of Eqs. (10) and (11) expressed in terms of <f> are 

Gi({& lifc })^ (6) =gj and G 2 {{b izk })^ b) < g 2 , (18) 

where 

Gi(K*}) = G 1 S({&, |Jfc })DA, 

G 2 ({^}) = G 2 S({6,, fc })DA, 
gi = gi - G lS ({6, )fc }) - GiS({&,, Jfc })Dp, 
g 2 = g 2 - G 2 s({6 i , Jfc }) - GiS({6 t(i fe})Dp. 

D.3 Optimization Over The Breakpoint Variables 

We have detailed how to determine *i%(xA^}A^l}) given a specification of all {{£>} and but 
how do we determine the best feasible choice of {&$} and {&$}? Before we answer that question, we have 
to examine the infeasibilities that might arise due to improper choices of {&• ^} and 

D.3.1 Elimination of {b i;k } Infeasibilities 

Most of the settings of or {b[*l} will not yield a feasible solution because Eq. (13) or Eq. (12) will be 

violated. The reason is that, for any fc, it is highly likely to get two segments x { k h) (b tth ) < xf ] < x[ b) (b l>k + 1) 
and x k b \b^ k ) < x[ b) < x k b) (b j>k + 1) that do not overlap (i.e., the optimizer will not be able to find a x {b) 
value that satisfies both bounds, and hence no feasible solution will be found). The same is also true for 
supplier response segment x^ . 

To eliminate those infeasibilities, we construct a super set of breakpoints {zu b k } such that 

4 6) = lK4 b) (^ = D,*«(6W = 2), • • • = 

i=l 

Hence each piecewise linear function f i>k will be described by a set of breakpoints zu k = {x fc (l), - - • , ar A («*.)} 
and function values f iyk (x k (l)), • ■ - J t , k {x k {K k )). These function values at the new breakpoints are lin- 
early interpolated from the old breakpoint-function pairs. After this reduction we have «* - 1 different 
b k e [1, K k - 1] values, which eliminates the need for constraint Eq. (13). 

The continuous supplier variables and constraint Eq. (12) are more problematic. Suppressing any 
dependence on we know 

4 s) = • • • , e^oei}) + £<m:m s) + e 4»4 6) - 

k<i k 
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If i^h and > * * ' y^Li} are known, then this expression allows us to bound the values x[ s) may take 
on, and thus bound the interval in which the feasible b [ "J lie. If we define the functions 

smax(>,y) if s > 0 \ smm(x,y) if s > 0 

smax(s,x i y) = < and smin($, x, y) = i 

jmin(x,y) if s < 0 |^max(a;,y) if s < 0 

then the minimum and maximum values that x[ s ^ may take on are: 

min(^>) = caftW,. ■ ■ , 6 W l} , + E^in(m«(^),x«( & W),^)(^ + l)) 

+ £ -in«(^), 4 6) + 1)) , (19) 

k 

and 

max(^) = ... ^ij,^}) +E sm -K s 2(< 5 2)>4 s) (^)>4 s) (< s i + D) 

+ E^WIdS).*?'©,*?'^! + 1)) (20) 

k 

respectively. The above equations imply that we can fully translate a complete specification of buyer variable 
breakpoint intervals to a corresponding set of feasible supplier variable breakpoint intervals. We will search 
over the breakpoints using a branch-and-bound algorithm, and will need to be able to infer the restriction on 
the values of feasible supplier breakpoint intervals given a partial specification of buyer breakpoint interval 
variables. This inference, which narrows the domains of the supplier breakpoint variables, is called con- 
straint propagation. Before describing how the branch-and-bound algorithm search over breakpoint interval 
variables is accomplished, we show how the constraint propagation may be obtained. 

Narrowing of Supplier Breakpoint Domains Through Propagation 

Let /3 ( k b) be a vector of integers specifying a set of breakpoint intervals for the Jfeth buyer variable. \0 k \ = 1 if 
the interval is uniquely determined, and \/3 k \ = K ( k b) - 1 if no restrictions on the interval have been specified. 
We index the Ith component of as / e [1, 4 &) - 1]. A set of breakpoints for all buyer dimensions 
is written as In analogous fashion, we define $\ 0$, I € [l,^ - 1], and {j3 k } for the supplier 

breakpoint interval variables. 

From a specification {^ 6) } we infer the feasible domains {p { k ] } using max(^ s) ) and mm(x[ 8) ) where: 

max(*W) = £ -ax {e^({6^}) + D)} + 

k b k e Pk J 

k<l b k E Pk J 



43 



and 

min(*i«>) = ^^^{^({6^}) + sm^rn^U^^^^ + 1))} + 

E J^K^}) + smin^^ftW),^^),^^^ + 1))} (22) 

k<l b k 

We begin with % = 1 (where no supplier variables contribute to x[ s ^) and determine the interval Ii = 
[min(4 s) ),ma^(2;[ s) )]. We then form the vector fi[ s) with components {p[ s J \[x[ s) ^ (/?i ; j + l)]nli # 0}. 

Equations (21) and (22) can then be iterated to find the interval J 2 = [min(a:^), max(x^)] using the 
previously calculated min(a^) and max(a;^). With I 2 known, the components of 0^ are determined as 
above. Iterating this procedure, we can find the potentially narrowed domains of all {0^}. 

We now return to the optimization problem: we need to search over all Yl k x EU ^ combinations of 
breakpoints to determine the combination that minimizes d\ (x opt (x, The simplest manner 

in which this accomplished is via exhaustive enumeration over all {b^} and {b^}. However, once all b^ 
have been specified, we need only search exhaustively over the corresponding {0k($)} inferred as above. 
There are, however, many ways to speed up the search over the buyer breakpoint intervals. We describe two 
methods. 

D.3.2 Local Search 

This method is extremely fast but not guaranteed to return the global optimum across all buyer breakpoint 
intervals. Here, initial random and feasible 19 {&j^}> and {b^} sets are first created. After the continuous 
optimization is carried out using the initial {6^}, and {6^}, a new {6^}, and {b^} set is generated from 
the old one by examining whether each generated x^\{b^}, {b^}) solution is on the boundaries of the 
interval x[ b) (b k ) < xf ] < x ( k b) (b k + 1) or x { k s \b k ) < x[ s) < x[ s) (b k + 1). If x[ b) is on the lower bound of the 
interval, bf^ = b[^ - 1; similarly, if it is on the upper bound, b k h) = b^ + 1. If is not on either end, 
then = b[ b K Of course this done while maintaining b^ in the interval 1 < b k b) < K k h) - 1. Supplier 
breakpoint variables are updated in a similar manner. The procedure is repeated until no decrease in distance 
is attained. 

D.3.3 Branch- an d-Bound Search 

Branch- and-bound is an efficient procedure used to exhaustively search a space and eliminate configurations 
that can be proven to be suboptimal. The theory behind branch-and-bound is well understood [5]. The 
idea is to structure the space as a tree of partially-completed possibilities and to explore the tree so that 
good complete solutions are constructed early on. Branches of the tree that can be shown to yield solutions 
19 Feasible in the sense that they do not violate the constraints in equations (10) , (11), (12), and (13). 
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which cannot be better than the currently best solution need not be explored. Heuristics indicate the order 
in which to explore the tree and provide a lower bound on the fitness of partially completed solutions for 
pruning of suboptimal branches.. 

In the present context, we carry out branch-and-bound search over the space of possible buyer breakpoint 
combinations. The tree is structured beginning from a node at which none of the b[ b) have been assigned. 
From this node, the tree branches into the *< 6) - 1 possible assignations for b\ b) . From any of these nodes, the 
tree branches into the possible values for bf\ The tree branches at each node to assign the remaining 
variables. At the leaves of the tree, all buyer breakpoint variables have been assigned. 

From any node, the search proceeds by choosing the most promising node to consider, as indicated by 
a variable ordering heuristic and considering the possible values for that variable in the order indicated 
by the value ordering heuristic. A complete solution and its accompanying distance is obtained when the 
search ends at a leaf of the tree. To avoid unnecessary exploration, a lower bound heuristic indicates when a 
partially completed solution (node) need not be expanded because all completions will have higher distance 
than the best located so far. We describe these three heuristics for search over buyer breakpoints. 
Lower Bound 

Given a set of buyer breakpoint interval vectors {^ 6) }, we can propagate to {p^} and the intervals 
[min(^ s) ),max(^ s) )]. A lower bound on the value d± can be obtained by ignoring the supplier capabilities 
and minimizing 

di(x) = (xW-/i)*C- 1 (xW- A i) (23) 

subject to 

Gix (s) - gi (24) 
G 2 x (s) < g 2 (25) 
mm{x[ s) ) < x^ < max(x[ 8) ), (26) 

The lower bound thus obtained is denoted by \b({p^}). 
Variable and Value Ordering 

The lower bound lb({/?( 6 )}) can be used to define an ordering over variables and the values any variable may 
take on. We define the functions L kyi which operate on the set {£ (6) } by locking the kth buyer breakpoint 
variable to p k . t . At any node we can determine the lower bounds lb(L fc) j{/?W}). We can select to explore 
the variables in order of increasing lb& where 
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i.e., the average over all possible completions for b[ j . 20 Within any given variable the (say k*) values can be 
explored in order of increasing lb(L k *j({p^})). If we are willing to do more computation, we can explore 
the branches in an interleaved order according to all the \b(L ki i({^})) . 

Pruning and Search 

During the branch-and-bound search, branches can be pruned if the node has lb(L k) i({0^})) greater than 
the distance of the lowest examined leaf node. 

Search begins at the root node where each pf> = [1, • • ■ - 1]. We applying the ordering heuristics 

above to move to a node where one breakpoint variable has been fixed. From that node, we move down the 
tree, fixing a value at a time until we reach a leaf. We then progressively backtrack up the tree, eliminating 
branches whose lower bounds are greater than the best leaf examined. 

Once the optimal set of buyer breakpoints has been determined, we search exhaustively over the supplier 
breakpoints feasible with this choice according to Equations (19) and (20). 

Improvements 

We have not yet exploited the fact that, from any node having {jSjj, 6 *}, the constraints in Equation (12) 
may narrow the domains of the buyer breakpoint variables deeper in the tree. The computational cost of 
narrowing the breakpoint domains may be warranted if the pruning thus obtained is sufficient. 

E Phase 2 Optimization 

In this section we detail an efficient algorithm to optimize Eq. (5): 

arg min di(x(x), >*) + d 2 {x) subject to C w (x) A C (s) (at). 

XT 

Our approach is a marriage of constraint programming and local search. 
E.l Local Search over ^-neighborhoods 

We begin with a definition. The cf-neighbors of a discrete configuration >c are all those configurations, V, 
which differ from >c in d attributes. Attribute % differs in two configurations if >q ^ ^. Thus the distance 
varies from 0 if the configurations are identical, up to n d if all attributes are different. The ^-neighborhood 
of h is the union of all d = 1 neighbors, d = 2 neighbors, up to d = D neighbors. The D-neighborhood 
around yc is indicated by Nd{k)- We note that, for any configuration, its n d neighborhood is the entire 
discrete space. 

Using the D-neighborhood, we apply a standard hillclimbing algorithm which respects all the discrete 
constraints. D is specified before the algorithm begins and serves as an efficient manner to trade off running 
20 We might also take the minimum over all possible completions \b k ({0^}) = min ; ^ \b(L ki i({0^}}) . 
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Table 5: Configurations examined in the D = 3 neighborhood search around (a) the initial configuration 
11111 and (b) the initial configuration 0+0+0 (duplicate search configurations are indicated in bold) . 



time with the quality of the final solution. The hillclimbing algorithm takes as input a starting configuration 
7t t and finds the configuration H t+1 e M D {>Ct) which minimizes di(x(xr),x) + d 2 (x) subject to C^fr) A 
C (s) (>*). The algorithm then iterates again starting from >f t+1 . The hillclimbing algorithm terminates when 
it finds a configuration h t which has a distance less than or equal to all of its D-neighbors. If D — n d , 
this is the global minimum. Given a function which exhaustively searches a D-neighborhood for the lowest 
distance feasible configuration, a D-neighborhood hillclimbing algorithm is easily implemented. 

There is one noteworthy point about the implementation of an efficient D-neighborhood hillclimber. 
Consider a D = 3 search over n d = 5 discrete variables. Without loss of generality, we label the initial 
configuration as 11111. In the D = 3 neighborhood search, the configurations examined are given in Table 
5(a) as indicated by the *s. Let us assume that the best neighbor in the D neighborhood of 11111 differs in 
components 2 and 4. We indicate this best neighbor as 1+1+1, which serves as the starting configuration 
for a new D — 3 search. The configurations examined in the second neighborhood search are given in 
Table 5(b). From Table 5(b) we see that the search branches 2, 7, and 10 are redundant, having already 
been performed in the D neighborhood search around 11111. 21 Thus these branches may be pruned in the 
second D neighborhood search. A search branch covers a previous search and may be pruned if the vector 
of positions of the wildcard characters covers the vector of positions at which the previous best solution was 
found. The binary relation covers is defined as 

x covers x >cl <^ yc V >c[ e *c* . 



x The duplication of search effort is most acute when D is large and the fittest neighbor is found very close by (small d). 
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The covers relation is transitive but not reflexive. Thus, to increase the efficiency, we can keep a list of the 
previous branches searched and ignored search branches that have already been covered. 

E.l.l Improvements 

Throughout the present discussion we have assumed that D remains constant through the hillclimbing 
process. There are good reasons why this should not be the case. With any reasonable landscape correlation 
structure, we know that lower distance neighbors are more likely to be found early in the hillclimbing process. 
Thus a reasonable heuristic is to start with a larger D early on and gradually lower it with time. 

E.2 Exhaustive Search of a ^-neighborhood 

To complete the discussion of the phase 2 optimization, we describe a computationally efficient way in which 
to find the lowest distance feasible configuration within a ^-neighborhood. We use a branch- and-bound 
method with heuristics tuned for the present problem. The factors tuned include the variable order in which 
we search, the order of the values any given variable may take, and a heuristic which says when certain 
branches of the search space may be safely ignored. 

We use heuristics which early on in the search generate low-distance configurations, and later prune parts 
of the search that we can prove must have higher distance than the best configuration generated so far. 

E.2.1 Heuristics And Lower Bounds 

We focus on obtaining good heuristics for the discrete contribution to the total distance, then later incorpo- 
rate the continuous and range contributions since these are simpler. 

At any point in the search, some variables have been assigned and some have not. We divide up the sum 
over the Z % ^{k z ,k^) into four contributions 

Q d Z w (x) = Q d ]p ]P w z w 3 Z i>J (>c l ,>€ J ) +Q d w t w 3 Z itj (x t ,x 3 ) 

% locked j locked % locked 3 floating 

i floating j locked % floating 3 floating 

= Z^i + Z\j + Z^\ + Zfj 

Since we are minimizing and all Z %y3 -{x u Xj), w u and Q d are positive, we can prune any search node which 
has Zi f i > than the distance of the lowest distance configuration found thus far, since additional terms will 
only add to the sum. 

We obtain a lower bound on the distance by lower-bounding each of the above four terms. Just taking 
the lowest value in the tables can yield differing values for any given variable, e.g., the values for a single 
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variable >c % \ 



k* = arg min Z ti3 (x t , x 3 ), = arg min Z 3i1 {k 3 ,h 1 ), and < = arg min Z M (x») 
may all be different. To correct this we write 

QdZw = Z ltl + Z f , { + ]r Z \}«) 

i floating 

where Zy\ and Zf )f are calculated as before and 

Z^H^^Qd^Zi^+Wi ]T w 3 {Z li3 {>c u ^) + z u {j< h ^))y 

j locked 

The best choice for the floating variable ^ is 

H ; = arg min z\ % }(>t % ). (27) 

In this way all the possible choices for >c { are at least consistent, and the true lower limit is still guaranteed 
to be higher than this (since the Z fjf contribution can only be higher). 

This approach also naturally suggests a variable ordering and a value ordering within any variable. We 
can try the variables in increasing order of their corresponding Z$(x*) so that variable i, which has the 
smallest Z$(xf), is tried first. For any given variable k % we can also order the possible value assignations 
for ^ according to the array of values Z$(x t ) so that the value K t = is tried first and the value which 
maximizes z\ % J{h z ) is tried last. 

With little extra work we can impose a consistency requirement raising the lower bound for Zf f . We 
decompose the sum as follows 

Zf, f = Qd WiWjZi^i,^) = Q d ^(^i(x t> ^) + Z h% (n 3 ,x % )). 

i,j floating: j^i 2jJ floating: j<i 

Instead of taking the minimum component- wise as previously suggested, i.e instead of 

Z tf = Qd ^2 ^ li ^[ W ^ W 3 Z t , 3 (H l ,>C 3 )], 
z,j floating 3 

we take the minimum so that 

z tf=Qd min^t^^,^) + Z 3 ^(x 3J Xi))]. (28) 

t,j floating: j<j 

which enforces pairwise consistency between variables i and j. 
The lower bound, lb, we have obtained is thus 

ib^^ + ^ + ^z/^o 
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where yt\ is determined from Eq. (27) and Z f * f is determined from Eq. (28). 
Integration Including Continuous and Range Contributions 

In this section we summarize the complete heuristic specification, including the effects of the continuous and 
range variables. As above, we use the definitions 

- Qd ^ w t WjZ it j(>a 9 >Cj) 

i y j locked 
j locked 

Z h = Qd rmn[w t Wj(Zi,j(xi,3€ 3 ) + Zj^x^Xj))] 

i y j floating. j <t 

so that a lower bound to the complete discrete contribution to the distance is Z\ t i 4- Zf f + J2 t Z* f (x*) where 
x* = arg min^ Z^(>a). 

We now extend this result to include the continuous and range variables. The continuous and range 
variables depend on the discrete variables through the mask. Since we expect there will be relatively few 
masks we exhaustively enumerate all the discrete possibilities. Consequently, we define d\ (boating flocked) 
as the distance contribution made by the continuous and range variables (the d\ contribution to Eq. (5)). 
We define the continuous + range contribution for any given floating variable i as 

= a m . m v d l (^floating flocked) 
noating\t 

where the min is taken over all floating variables other than i and where floating variable i is locked to value 

Since each continuous/range/discrete distance contributes additively to the total distance, we can rank 
the possible values for any given floating variable % according to Z$(x % ) + d^(v % ). The best (according to 
this heuristic) setting for floating variable % is 

= arg min + df {v % )\ 

The heuristic then suggests a variable ordering in order of increasing z\ % }{y<l) + d[ l \v^) so that we first 
select the floating variable with the lowest sum and last select the variable with the highest sum. 

Finally, we turn to the lower bound on the distance. Noting that the continuous and range distance 
contribution does not have the form of a sum over the dimensions (as the discrete contribution does) , we 
write the lower bound, lb, as 22 

lb = Z ltl + ^(2| ( }K) + 6 ltt . + Z* u 

z 

where z* = arg min, (z{*)(*c*) + d[ l \v*)) (i.e., i* labels the floating dimension with the least Z$ + d[ 1 ^ 
contribution to the distance). 

22 $x t x' is tne Kronecker delta function which is 1 if and only if x ~ x } . 
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E.2.2 Improvements 

We point out that an extensive literature exists on parallel branch and bound algorithms. Many of the 
techniques developed are applicable to the present problem. 
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